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Abstract. Our contribution is two-folded. First, starting from the known fact that every real 
skew-Hamiltonian matrix has a real Hamiltonian square root, we give a complete characterization of 
the square roots of a real skew-Hamiltonian matrix W. Second, we propose a structure-exploiting 
method for computing square roots of W. Compared to the standard real Schur method, which 
ignores the structure, our method requires significantly less arithmetic. 
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1. Introduction. Given A 6 C" xn , a matrix X for which X 2 = A is called a 
square root of A. The matrix square root is a useful theoretical and computational 
tool, one of the most commonly occurring matrix functions. See jTTTl H51 IT51 UM H71 12U] . 

The theory behind the existence of matrix square roots is nontrivial and the 
feature which complicates this theory is that in general not all the square roots of a 
matrix A are functions of A. See [51 120]. 

It is well known that certain matrix structures can be inherited by the square 
root. For example, a symmetric positive (semi) definite matrix has a unique symmetric 
positive (semi) definite square root |19) . The square roots of a centrosymmetric matrix 
are also centrosymmetric |23j . A nonsingular M-matrix has exactly one M-matrix 
as a square root. For an iJ-matrix with positive diagonal elements there exists one 
and only one square root which is also an iJ-matrix with positive diagonal elements 
[21) . The principal square root of a centrosymmetric H- matrix with positive diagonal 
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elements is a unique centrosymmetric if -matrix with positive diagonal entries [22j. 
Any real skew-Hamiltonian matrix has a real Hamiltonian square root [8 . In this 
paper we characterize such square roots. 

For general matrices, an attractive method which uses the Schur decomposition is 
described by Bjorck and Hammarling [4 but may require complex arithmetic. Higham 
[12] presented a modification of this method which enables real arithmetic to be used 
throughout when computing a real square root of a real matrix. This method has 
been extended to compute matrix pth roots [27j and general matrix functions [7] . 

It is a basic tenet in numerical analysis that structure should be exploited allow- 
ing, in general, the development of faster and/or more accurate algorithms [3 [25] . 
We propose a structure-exploiting method for computing square roots of a real skew- 
Hamiltonian matrix W which uses the real skew-Hamiltonian Schur decomposition 
and requires significantly less arithmetic. 

We give some basic definitions and establish notation in Section 2. A description 
of the real Schur method and some results concerning the existence of real square 
roots are presented in Section 3. In Section U we characterize the square roots of 
a nonsingular W in a manner which makes clear the distinction between the square 
roots which are functions of W and those which are not. In Section [S] we present our 
algorithms for the computation of skew-Hamiltonian and Hamiltonian square roots. 
In Section [6] we give results of numerical experiments. 



2. Definitions and preliminaries results. 



2.1. Square roots of a nonsingular matrix. Given a scalar function / and a 
matrix A £ C nx ™ there are many different ways to define f(A), a matrix of the same 
dimension of A, providing a useful generalization of a function of a scalar variable. 



It is a standard result that any matrix A S 
canonical form 



can be expressed in the Jordan 



Z^AZ = J = diag( Ji, J 2 , . . . , J 3 



Jk — Jk{^k) 



where Z is nonsingular and mi + to 2 
up to the ordering of the blocks Ji. 



A* 



1 



1 

Afe. 



pj> 



-im fc xm k 



(2.1) 
(2.2) 



n. The Jordan matrix J is unique 



Denote by Ai , . . . , A s the s distinct eigenvalues of A and let rii be the order of the 
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largest Jordan block in which A.; appears. The function / is said to be defined on the 
spectrum of A if the values 

/ (j) (Ai)> 3 = 0,...,ni - 1, i=l,...,s, 
exist. These arc called the values of the function f on the spectrum of A. 

The following definition of matrix function defines f(A) to be a polynomial in the 
matrix A completely determined by the values of / on the spectrum of A. See [HI p. 
407 ff.]. 

Definition 2.1 (matrix function via Hermite interpolation). Let / be defined 
on the spectrum of A £ C nxn . Then 

f(A) :=p(A) 

where p is the polynomial of degree less than J^i n i which satisfies the interpolation 
conditions 



P U) (Xi) = f U) (Xi), .7=0, 



1, i = l,...,s. 



There is a unique such p and it is known as the Hermite interpolating polynomial. 

Of particular interest here is the function f(z) = z 1 / 2 which is certainly defined 
on the spectrum of A if A is nonsingular. However, the square root function of A, 
/(A), is not uniquely defined until one specifies which branch of the square root is 
to be taken in the neighborhood of each eigenvalue A,; . Indeed, Definition 12.11 yields 
a total of 2 s matrices f(A) when all combinations of branches for the square roots 
f(Xi), i = 1, . . . , s, are taken. It is natural to ask whether these matrices are in fact 
square roots of A, that is, do we have f(A)f(A) = A? Indeed, these matrices, which 
are polynomials in A by definition, are square roots of A. See [HI [20]. However, these 
square roots are not necessarily all the square roots of A. 



To classify all the square roots of a nonsingular matrix A £ C" 
following result concerning the square roots of a Jordan block. 



we need the 



Lemma 2.2. For A& ^ the Jordan block Jfc(Afc) in (|2 .2[) has precisely two upper 
triangular square roots 



L^=L ( i\\ k ) = 



/(A fc ) /'(A fc ) 

/(A*) /'(A*) 



/(A fc ) 



/ ( '" fc ' 1) (A fc ) 

(m k -l)\ 
/ (mfc - 2) (A fc ) 
(m k -2)\ 



/'(Afc) 
/(Afc) 



J = 1,2, (2.3) 
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where /(A) = A 1 / 2 and the superscript j denotes the branch of the square root in the 
neighborhood of Xk ■ Both square roots are functions of Jk ■ 

We will restrict our attention to matrices with real entries and to investigate the 
real square roots of a real matrix we need to understand the structure of a general 
complex square root. The following results allow us to obtain a useful characterisation 
of the square roots of a nonsingular matrix A £ C" x ™ which are functions of A. See 

Theorem 2.3. Let A £ C»xn jj e nonsingular and have the Jordan canonical 
form (|2.2p . Then all square roots X of A are given by 

X = ZU diag (L^\L^\. . . , L^) U~ l Z-\ 

where jk — 1 or jk — 2 and U is an arbitrary nonsingular matrix which commutes 
with J. 

The following result extends Theorem 12.31 

Theorem 2.4. Let the nonsingular matrix A £ C nxn have the Jordan canonical 
form (|2.2p and let s < p be the number of distinct eigenvalues of A. Then A has 
precisely 2 s square roots which are functions of A, given by 

X 3 =Z^g(L^\L^\...,L^)z~\ j = l,...,2', (2.4) 

corresponding to all possible choices of j\, . . . , j p , jk — 1 or jk = 2, subject to the 
constraint that ji = jk whenever A^ = A^ . 

If s < p, A has square roots which are not functions of A; they form parametrized 
families 

Xj(U) = ZUdmg (L[ jl) ,L^\ . . U- l Z-\ j = 2 s + 1, . . . , 2", (2.5) 

where jk — 1 or jk — 2, U is an arbitrary nonsingular matrix which commutes with J , 
and for each j there exist i and k, depending on j , such that Xi — Xk while ji ^ jk . 

Proofs of these theorems and a description of the structure of the matrix U can 
be found in [10] . Note that formula in (|2.4I) follows from the fact that all square roots 
of A which are functions of A have the form 

f(A) = fiZJZ- 1 ) = Zf^Z- 1 = Zdiag (f(J k ))Z-\ 

and from Lemma 12.21 The constrain on the branches {ji} follow from Definition 12.11 
The remaining square roots of A (if any), which cannot be functions of A, are given 
by QMS). 
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Theorem 12.41 shows that the square roots of A which are functions of A are 
"isolated" square roots, characterized by the fact that the sum of any two of their 
eigenvalues is nonzero. On the other hand, the square roots which are not functions 
of A form a finite number of parametrized families of matrices: each family contains 
infinitely many square roots which share the same spectrum. 



Some interesting corollaries follow directly from Theorem 12.41 

Corollary 2.5. If Xk ^ 0, the two square roots of Jk(Xk) given in Lemma 12.21 
are the only square roots of Jk (Xk ) ■ 

Corollary 2.6. If A £ C" xn j s nonsingular and in its Jordan canonical form 
(|2.2[) each eigenvalue appears in only one Jordan block, then A has precisely 2 P square 
roots, each of which is a function of A. 

The final corollary is well known. 

Corollary 2.7. Every Hermitian positive definite matrix has a unique Hermi- 
tian positive definite square root. 

2.2. Hamiltonian and skew-Hamiltonian matrices. Hamiltonian and skew- 
Hamiltonian matrices have properties that follow directly from the definition. 



Definition 2.8. Let J = 



/ 

-/ 



, where / is the identity matrix of order n. 



(1) A matrix H e R 2nx2n is said to be Hamiltonian if HJ = (HJ) T . 
Equivalently, H can be partitioned as 



H 



A G 

F -A T 



G = G T , F — F T , A,G,F ei" x ". (2.6) 



(2) A matrix W £ R 2 " x2n is said to be skew-Hamiltonian if WJ = -{WJ) T . 
Likewise, W can be partitioned as 



W 



A G 
F A T 



G = —G T , F = —F T , A,G,FeR nxn . (2.7) 



These matrix structures induce particular spectral properties for H and W . No- 
tably, the eigenvalues of H are symmetric with respect to the imaginary axis and the 
eigenvalues of W have even algebraic and geometric multiplicities. 

Definition 2.9. 

(1) A matrix S G JR 2nx2n is said to be symplectic if SJS T = J. 

(2) A matrix U € j]j2nx2n j g ^ Q ^ e or ihogonal- symplectic if U is 
orthogonal and symplectic. Any matrix belonging to this group can be par- 
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titioned as 



U = 



Ut u 2 

-U 2 Ui 



where Ui G 



i = 1,2. 



Hamiltonian and skew-Hamiltonian structures are preserved if symplectic sim- 
ilarity transformations are used; if H is Hamiltonian (skew-Hamiltonian) and S is 
symplectic, then S~ 1 HS is also Hamiltonian (skew-Hamiltonian). In the interest of 
numerical stability the similarities should be orthogonal as well. 

The first simplifying reduction of a skew-Hamiltonian matrix was introduced by 
Van Loan in [24] . But first we recall the real Schur decomposition pj] . 



Theorem 2.10 (real Schur form). If A 6 
onal matrix Q such that 



1 , then there exists a real orthog- 



Q AQ = R 



Rn 



R\2 

R22 



Rim 
R2m 

Rm/tn. 



(2.8) 



where each block Ru is either lxl or 2 x 2 with complex conjugate eigenvalues A.; 
and Xi, Xi ^ Xi (R is in quasi-upper triangular form). 

In [51] it was shown that any skew-Hamiltonian W can be brought to block-upper- 
triangular form by an orthogonal-symplectic similarity. Actually, we can explicitly 
compute an orthogonal-symplectic matrix U such that 



U T WU 



Wi W 2 

o 



where W2 



-W2 and Wi is upper Hessenberg (a matrix is upper Hessenberg if all 
entries below its first subdiagonal are zero). This is called the symplectic Paige/Van 
Loan (PVL) form. 

Subsequently, if the standard QR algorithm is applied to W% producing an or- 
thogonal matrix Q and a matrix in real Schur form N± so that 

Wi = QNiQ T , 

we attain the rea] skew-Hamiltonian Schur decomposition of W, 



U T WU = 



O N[ 



(2.9) 
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where U = U 



Q 
O Q 



and N 2 = Q T W 2 Q. 



Lemma 2.11 (real skew-Hamiltonian Schur form). Let W e R 2nx2 " be skew- 
Hamiltonian. Then there exists an orthogonal matrix 



U 



U x U 2 
-U 2 Ux 



such that 



U T WU 

and N\ is in real Schur form. 



'N t N 2 
N[ 



No 1 



-No 



(2.10) 



In |26[ Theorem 5.1] we can find a result concerning the real Hamiltonian Schur 
decomposition. 

Theorem 2.12 (real Hamiltonian Schur form). Let H e K 2 ?ix2n ^ e Hamiltonian. 
If H has no nonzero purely imaginary eigenvalues, then there exists an orthogonal 
matrix 



U 



Ux U 2 
-Ui Ux 



Ux,U 2 g 



such that 



U T HU 

and Hx is in real Schur form. 



Hx 




H 2 
-Hi 



Ho 



(2.11) 



In this article we are interested in the computation of a real square root of a 
real skew-Hamiltonian matrix W and to discuss square roots of a skew-Hamiltonian 
matrix we need to consider a variant of the Jordan canonical form (I2.2[) when the 
matrix A is real. In this case, all the nonreal eigenvalues must occur in conjugate 
pairs and all the Jordan blocks of all sizes (not just lxl blocks) corresponding to 
nonreal eigenvalues occur in conjugate pairs of equal size. 

For example, if A is a nonreal eigenvalue of the real matrix A, and if J2(A) appears 
in the Jordan canonical form of A with a certain multiplicity, J 2 (X) must also appear 
with the same multiplicity. See JTH1 p. 150 ff.]. The block matrix 



■MA) 

o 







"A 


1 
















A 








J 2 (A)_ 










A 


1 






.0 








A 



(2.12) 
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is permutation-similar (interchange rows and columns 2 and 3) to the block matrix 



A 





1 


0" 











A 





1 




D(X) 


/ 








A 










D(X)_ 











A. 









D(X) 



A 
O A 



(2.13) 



Each block D(X) is similar to the matrix 



SD(X)S~ 1 



a b 

—b a 



C(a,b), S 



1 -1 



(2.14) 



where A, A = a ± ib, a, b G KL b / 0. Thus, every block pair of conjugate 2x2 Jordan 
blocks (|2.12[) with nonreal eigenvalue A is similar to a real 4x4 block of the form 



a 


b | 


l 


0" 








-b 


a i 





1 




C(a,b) 


I 





i 


a 


b 




O 


C(a,b) 





o ! 


-b 


a . 









In general, every block pair of conjugate k x k Jordan blocks with nonreal A, 

J fc (A) O 
O J fc (A) 

is similar to a real 2k x 2k block matrix of the form 



(2.15) 



C(a,b) I 

C(a,b) I 



C(a,b) 



I 

C(a,b) 



(2.16) 



We call Ck(a,b) a real Jordan block. These observations lead us to the real Jordan 
canonical form. 

Theorem 2.13. Each matrix A e M ,ixn is similar (via a real similarity trans- 
formation) to a block diagonal real matrix of the form 

Cm ) 



Jr 



Cn v bp) 



Jn p+1 (Xp+l) 



(2.17) 
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where Xk = dk + ibk, ifc, l)fc £l, k = 1, . . . ,p, is a nonreal eigenvalue of A and Xk, 
k = p + 1, . . . ,p + q, is a real eigenvalue of A. Each real Jordan block C nk (a^, bk) 
is of the form (|2.16j) and corresponds to a pair of conjugate Jordan blocks J nk {Xk) 
and J nk (Xk) for a nonreal Xk in the Jordan canonical form of A (|2.1[) . The real 
Jordan blocks J nk (Xk) are exactly the Jordan blocks in (|2.ip with real Xk- Notice that 
2(m + • • • + n p ) + (n p +i + • • • + ?ip+<j) = n. We call Jr a real Jordan matrix of order 
n, a direct sum of real Jordan blocks. 

In [8] it is shown that every real skew-Hamiltonian matrix can also be reduced to 
a real skew-Hamiltonian Jordan form via a symplectic similarity. See also [S] . 

Lemma 2.14. £H Theorem 1] For every real skew-Hamiltonian matrix W G 
K 2nx2 ™ there exists a symplectic matrix ^ e R 2 " x211 such that 



where Jr £ M nx " is in real Jordan form (| 2 . 1 T[) and is unique up to a permutation of 
real Jordan blocks. 

In principle, a real square root of W can be obtained by the general method 
devised by Higham [T^]. Such method, however, does not exploit the structure of W. 
The method we propose exploit the skew-Hamiltonian structure of W and it also uses 
the real Schur method. To make our presentation simpler we decided to present the 
details of the real Schur method. See [12j p. 412 ff.]. 

3. An algorithm for computing real square roots. 

3.1. The Schur method. Bjorck and Hammarling [4 presented an excellent 
method for computing a square root of a matrix A. Their method first computes a 
Schur decomposition 



where Q is unitary and T is upper triangular [llj . and then determines an upper 
triangular square root Y of T with the aid of a fast recursion. A square root of A is 
given by 



A disadvantage of this Schur method is that if A is real and has nonreal eigenvalues, 
the method needs complex arithmetic even if the square root which is computed should 
be real. When computing a real square root it is obviously desirable to work with 
real arithmetic; depending on the relative costs of real and complex arithmetic on a 



Jr 



(2.18) 




Q*AQ = T 



X = QYQ*. 
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given computer system, substantial computational savings may occur, and moreover, 
a computed real square root is guaranteed. 

Higham described a generalization of the Schur method which enables real arith- 
metic to be used throughout when computing a real square root of a real matrix. In 
Section [3731 we present this method. First we give, for a nonsingular matrix A G M nx ", 
conditions for the existence of a real square root, and for the existence of a real square 
root which is a polynomial in A. 

3.2. Existence of real square roots. The following result concerns the ex- 
istence of general real square roots - those which are not necessarily functions of 
A. 

Theorem 3.1. Let A G R nx ™ be nonsingular. A has a real square root if and 
only if each elementary divisor of A corresponding to a real negative eigenvalue occurs 
an even number of times. 

Theorem 13.11 is mainly of theoretical interest, since the proof is nonconstructive 
and the condition for the existence of a real square root is not easily checked compu- 
tationally We now focus attention on the real square roots of A G M nx ™. The key to 
analysing the existence of square roots of this type is the real Schur decomposition. 

Suppose that A G M rtxrl and that / is defined on the spectrum of A and consider 
the real Shcur form of A in (|2.8j) . Since A and R in (|2.8|) are similar, we have 

f(A) = Qf(R)Q T , 

so that f(A) is real if and only if 

Z = f(R) 

is real. It is not difficult to show that Z inherits R's quasi-upper triangular structure 
and that 

Zu = f(Rii), i = l,...,m. 

If A is nonsingular and / is the square root function, then we have 

Z 2 = R and X 2 = A with X = QZQ T . 

The whole of Z is uniquely determined by its diagonal blocks. To see this equate 
blocks in the equation Z 2 = R to obtain 



j 

^ ZikZkj = Rij, j > i- 

k—i 
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These equations can be recast in the form 

Zf i = R i i, i=l,...,m (3.1) 

3-1 

ZuZij + ZijZjj = Rij — S ' ZikZkj, j = i + 1, . . . , to. (3-2) 

fc=i+l 

Thus, if the diagonal blocks Za are known, (|3.ip provides an algorithm for computing 
the remaining blocks Zij of Z along one superdiagonal at a time in the order specified 
by j — i = 1,2,. . . ,m — 1. The condition for the Sylvester equation (I3.2[) to have a 
unique solution Zij is that Zjj and — Zjj have no eigenvalue in common |ll[|20j . This 
is guaranteed because the eigenvalues of Z are fik = f(Xk) and for the square root 
function f(Xi) — —f(Xj) implies that Xi = Xj = 0, contradicting the nonsingularity 
of A. 

From this algorithm for constructing Z from its diagonal blocks we conclude that 
Z is real and hence f(A) is real, if and only if each of the blocks Za = f(Ru) is real. 

We now examine the square roots f(A) of a 2 x 2 matrix A with complex conjugate 
eigenvalues. Since A has 2 distinct eigenvalues, it follows from 

Corollary 12.61 that A has four square roots which are all functions of A. Next lemma 
says about the form os these square roots. 

Lemma 3.2. Let A £ R 2x2 have complex eigenvalues X,X ~ 9 ± where /x ^ 0. 
Then A has four square roots, each of which is a function of A. Two of the square 
roots are real, with complex conjugate eigenvalues, and two are pure imaginary, having 
eigenvalues which are not complex conjugate. 

Using this lemma it can be proved 

Theorem 3.3. Let A £ W nxn be nonsingular. If A has a real negative eigenvalue, 
then A has no real square roots which are functions of A. 

If A has no real negative eigenvalues, then there are precisely 2 r+c real square 
roots of A which are functions of A, where r is the number of distinct real eigenvalues 
of A and c is the number of distinct complex conjugate eigenvalue pairs. 

For proofs of Lemma 13.21 and Theorem 13.31 see [l2l p. 414 ff.]. 

It is clear from Theorem 13.11 that A may have real negative eigenvalues and yet 
still have a real square root; however, as Theorem 13.31 shows . the square root will not 
be a function of A. 

3.3. The real Schur method. The ideas of the last section lead to a natural 
extension of Bjorck and Hammarling's Schur method for computing in real arithmetic 
a real square root of a nonsingular A 6 R nx ™. This real Schur method begins by 
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computing a real Schur decomposition (|2.8j) . then computes a square root Z of R from 
equations (|3.1|) and (|3.2|) , and finally obtains a square root of A via the transformation 
X = QZQ T . 

The solution of Equation (I3.1[) can be computed efficiently in a way suggested by 
the proof of Lemma [3.21 12, p. 417]. The first step is to compute 9 and /i, where 
A = 6 + ifi is an eigenvalue of the matrix 



T21 r 2 2 



Next, a and /3 such that (a + i/3) 2 = 9 + i/j, are required. Finally, the real square roots 
of Ru are given by 

Za = ±(aI+±- (Ru -91) J. (3.3) 
If Zu is of order p and Zjj is of order q, Equation (|3 . 2[) can be written as 

(7 9 <g> Z lt + Zj } <g> 7 p ) col(Zy) = col f% - J2 Z 'kZ k j] , J > i (3.4) 

V fc=i+i / 

where eg) is the Kronecker product and col(M) denotes the column vector formed by 
taking columns of M and stacking them atop one another from left to right. The 
linear system (|3.4[) is of order pq — 1, 2 or 4 and may be solved by standard methods. 

Note that to conform with the definition of f(A) we have to choose the signs in 
(|3.3j) so that Zu and Zjj have the same eigenvalues whenever Ru and Rjj do; this 
choice ensures simultaneously the nonsingularity of the linear systems (|3.4|) . 

Any of the real square roots f(A) of A can be computed in the above fashion by 
the real Schur method. 

Algorithm 1 [Real Schur method] 

1. compute a real Schur decomposition of A, 

A = Q T RQ- 

2. compute a square root Z of R solving the equation Z 2 = R via 

Z 2 = R u , 1 < % < m, 

(l q ® Z u + Zjj ® I p ) col(%) = col ( Rij - V Z lk Z k3 V j > i 



[block fast recursion] 
3. obtain a square root of A, X — QZQ T . 



Square roots of real skew-Hamiltonian matrices 



13 



The cost of the real Schur method, measured in floating point operations (flops) 
may be broken down as follows. The real Schur factorization costs about 15n 3 flops 
[TTj . The computation of Z requires n 3 /6 flops and the formation of X = QZQ T 
requires 3n 3 /2 flops [HI p. 418]. Only a fraction of the overall time is spent in 
computing the square root Z. 

4. Square roots of a skew-Hamiltonian matrix. In this section we present 
a detailed classification of the square roots of a skew-Hamiltonian matrix W G K 2 nx2n 
based on its real skew-Hamiltonian Jordan form (|2.18p . 



Jr 



(4.1) 



where Jr 6 R" xn is in real Jordan form (|2.17|) and ^ is a sympletic matrix. First, 
we will discuss the square roots of Jr. 

According to Lemma [2~2l for ^ a canonical Jordan block Jk(^k) has precisely 
two upper triangular square roots given by (|2.3p . As a corollary it follows that 



Corollary 4.1. For a real eigenvalue A& ^ 0, the Jordan block Jfc(A^) in (12.21) 
has precisely two upper triangular square roots which are real, if A > 0, and pure 
imaginary, if X < 0. Both square roots are functions of </&. 

To fully characterize the square roots of a real Jordan block CV-(a, 6), we first 
examine the square roots of a 2 x 2 block C(a, b) corresponding to the nonreal eigen- 
values A, A = a ± ib, a, b 6 R. According to Lemma l3~2l the real Jordan block C(a, 6) 
in (12. 14)) has four square roots, each of which is a function of C(a,b). Two of the 
square roots are real and two are pure imaginary. 



Lemma 4.2. The real Jordan block Ck{a,b) in (I2.16P has precisely four block 
upper triangular square roots 



F, 



F Fx ■ 
F Fx 



F 



F k -i 

Fk-2 

F 



(4.2) 



where F is a square root of C(a, b) and Fi, i = 1, . . . , k — 1, are the unique solutions 
of certain Sylvester equations. The superscript j denotes one of the four square roots 
of C(a,b). These four square roots Fl o) are functions of Ck(a,b), two of them are 
real and two are pure imaginary. 



Proof. Since Ck(a,b) has 2 distinct eigenvalues and the Jordan form (|2.15|) has 
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p = 2 blocks, from Corollary 12.61 we know that Cfc(a, b) has four square roots which 
are all functions of A. 

Let X be a square root of Cfc(a, b) (k > 1). It is not difficult to sec that X inherits 
Cfe(a, b) block upper triangular structure, 



X = 



X22 



Xkk. 



where Xij are all 2 x 2 matrices. Equating blocks in the equation 



we obtain 



X 2 = C k (a,b) 



X% = C{a,b), i = l,...,k, 
XnXij+i + Xi^ + iXi + \^ + i = I2, i = 1, . . . , k — 1 

3-1 

XuXij + XijXjj = — XuXij, j = i + 2, . . . , k. 

l—i+l 



(4.3) 



(4.4) 
(4.5) 

(4.6) 



The whole of X is uniquely determined by its diagonal blocks. If F is one square root 
of C(a,b), from (|4.4I) and to conform with the definition of f(A) (the eigenvalues of 
Xa must be the same), we have 



X U = F, i = l,.. 



(4.7) 



Equations (|4.5[) and (|4.6[) are Sylvester equations and the condition for them to 
have a unique solution Xij is that Xu and —Xjj have no eigenvalues in common and 
this is guaranteed. 

From (14.5[) we obtain the blocks X^ along the first superdiagonal and (|4.7[) forces 
them to be all equal, say Fi, 

X\2 = X2Z = . . . = Xh-i t k = Ft- 

This implies that the other superdiagonals obtained from (|4.6p are also constant, say 
Fj-i, j = 3,...,fc, 

X%j = X2J+1 = . . . = Xk-j+i,k = Fj-x, j = 3, . . . , fc. 



Thus, since there are only exactly four distint square roots of C(a,b) which are 
functions of C(a, 6), F = F"\ I = 1,...,4, it follows that Ck{a,b) will also have 
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precisely four square roots which are functions of Ck(a,b). If F is real then Fj—i, 
j = 2, . . . , fc, will also be real. If F is pure imaginary it can also be seen that i^-i, 
j = 2, . . . , k, will be pure imaginary too. □ 

Next theorem combines Corollary 14.11 and Lemma 14.21 to characterize the square 
roots of a real Jordan matrix Jr. 



Theorem 4.3. Assume that a nonsingular real Jordan matrix Jr in (|2.17|) has p 
real Jordan blocks corresponding to c distinct complex conjugate eigenvalue pairs and 
q canonical Jordan blocks corresponding to r distinct real eigenvalues. 



Then Jr has precisely 2 2c+r square roots which are functions of Jr, given by 



Xj = diag 



ptip) r(n) 
1 ii P ' iip+i ' 



i = i, - - • , 2 



2c+r 



(4.8) 



corresponding to all possible choices of ji, . . . , j p , jk = 1,2,3 or 4- cind i±,...,i q , 
ik = 1 or 2, subject to the constraint that ji — jk and i; = ik whenever A; = Xk- 

If c + r < p + q, then Jr has square roots which are not functions of Jr and they 
form 2 2p+q — 2 2c+r parameterized families given by 



J =2 



p+i ' 

_ 2c+r 



i, 



(4.9) 



2 2 P+9 



where jk = 1,2,3 or 4 and ik — 1 or 2, Q is an arbitrary nonsingular matrix which 
commutes with Jr and for each j there exist I and k depending on j , such that A/ = Xk 
while ji ^ j k or i x ^ i k . 

Proof. The number of distinct eigenvalues is s = 2c+r and, according to Theorem 
12.41 Jr has precisely 2 s = 2 2c+r square roots which are functions of Jr. All square 
roots of Jr which are functions of Jr satisfy 



f(c ni ) 



f(Cn p ) 



f(Jn p+1 , 



f{Jn p+q ) 



and, according to Lemma [2721 and Lemma 14721 these are given by (|4.8[) . The constraint 
on the branches {jk} and {ik} comes from Definition 12.11 of matrix 
function. The remaining square roots of Jr, if they exist, cannot be functions of 
Jr. Equation (|4.9p derives from the second part of Theorem [231 D 
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From Theorem 14. 3[ Lemma 14.21 and Corollary 14.11 the next result concerning the 
square roots of A which are functions of A follows immediately. 

Corollary 4.4. Under the assumptions of Theorem 14. 3\ 

(1) if Jr has a real negative eigenvalue, then Jr has no real square roots which 
are functions of Jr; 

(2) if Jr has no real negative eigenvalues, then Jr has precisely 2 c+r real square 
roots which are functions of Jr, given by (|4.8j) with the choices of ji, . . . , j p 
corresponding to real square roots Fni\ ■ ■ ■ , Fn p p ^ ; 

(3) if Jr has no real positive eigenvalues, then Jr has precisely 2 c+r pure 
imaginary square roots which are functions of Jr, given by (|4.8j) with the 
choices of ji , . . . , j p corresponding to pure imaginary square roots 
pO'i) pip 

r n% )***3 r n v ■ 

Now we want to use all these results to characterize the square roots of a real 
skew-Hamiltonian matrix W . 

Theorem 4.5. Let W G R 2 ™ x2 ™ he a nonsingular skew-Hamiltonian matrix with 
the real skew-Hamiltonian Jordan form in (|4.ip . Assume that Jr has p real Jordan 
blocks corresponding to c distinct complex conjugate eigenvalue pairs and q canonical 
Jordan blocks corresponding to r distinct real eigenvalues. 

Then W has precisely 2 2c+r square roots which are functions of W , given by 

Yj =*diag(X,,X, T )*-\ j = l,...,2 2c+r , (4.10) 

where Xj is a square root of Jr given in |^.($[ ). 

W has always square roots which are not functions of W and they form A 2p+q — 
2 2c+r parameterized families given by 

YjiQ) =t>edi&g(x j ,Xf^Q- 1 t>- 1 , j = 2 2c+r + l,...,4 2p+q , (4.11) 

where 

Xi = diag (f£ ),..., Fg) , , ... , IM) , 

Xj = diag (F^\ Fg»\L%£\. . . , L<$ q ) , 

jk = 1,2,3 or 4 omd ik = 1 or 2, is an arbitrary nonsingular matrix which commutes 
with diag (Jr, JjQ and for each j there exist I and k depending on j, such that A; = 
while ji ^ j k or % x ^ i k . 

Notice that W has s = 2c + r distinct eigenvalues corresponding to 2(p + q) real 
Jordan blocks and 2(2p + q) canonical Jordan blocks. We always have s < 2p + q and 
so s < 2(2p + q). Thus, there are always square roots which are not functions of W. 
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Proof. This result is a direct consequence of Theorem 12.41 and Theorem 14.31 
Notice that if Xj in (|4.8[) is a square root of Jr then diag [Xj , Xj) is a square root 
of diag (J R , JjQ. Thus, 



W = # 



Jr 



XJ. 



$-1 



Thus, 



X; 



is a square root of W . Since Xj is a function of Jr, Yj is a function of diag (Jr, Jr). 
This proves the first part of the theorem. 



The second part follows from the second part of Theorem 12.41 and the fact that 
diag (Xj,XjJ is a square root of diag (Jr, Jr). □ 

It is easy to verify that if W is a skew-Hamiltonian matrix, then W 2 is also 
skew-Hamiltonian. This implies that any function of W, which is a polynomial by 
definition, is a skew-Hamiltonian matrix. Thus, all the square roots of W which are 
functions of W are skew-Hamiltonian matrices. The following result refers to the 
existence of real square roots of a skew-Hamiltonian matrix. 

Corollary 4.6. Under the assumptions of Theorem 14.51 the following state- 
ments hold: 

(1) if W has a real negative eigenvalue, then W has no real skew-Hamiltonian 
square roots; 

(2) if W has no real negative eigenvalues, then W has precisely 2 c+r real skew- 
Hamiltonian square roots which are functions of W , given by (|4.10|) with the 
choices of ji , . . . , j p corresponding to real square roots F { n { l] ,...,F { n i p) . 

It is clear from Theorem 13.11 that W may have real negative eigenvalues and yet 
still have a real square root; however, the square root will not be a function of W . 

In (8j Theorem 2] it is shown that 

Lemma 4.7. Every real skew-Hamiltonian matrix W has a real Hamiltonian 
square root. 
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The proof is constructive and the key step is based in Lemma 12.141 - we can 
bring W into a real skew-Hamiltonian Jordan form (|4.1j) via a sympletic similarity. 
Further, it is shown that every skew-Hamiltonian matrix W has infinitely many real 
Hamiltonian square roots. 

The following theorem gives the structure of those real Hamiltonian square roots. 

Theorem 4.8. Let W £ l 2liX2,rl be a nonsingular skew-Hamiltonian matrix and 
assume the conditions in Theorem 14.51 



(1) If W has no real negative eigenvalues, then W has real Hamiltonian square 
roots which are not functions of W and they form 2 p+q parameterized families 
given by 



Yj(Q) = *6 diag (Xj, -X S T ) e- 1 *- 1 , j = l,..., 2 p+q , 



(4.12) 



where Xj denotes a real square root of Jr and is an arbitrary nonsingular 
symplectic matrix which commutes with diag (Jr, Jr)- 
(2) IfW has some real negative eigenvalues, then W has real Hamiltonian square 
roots which are not functions of W and they form 2 p+q parameterized families 
given by 



Yj(Q) = m 



Xj 
K, 



K 



3 

Xj 



(4.13) 



where Xj is a square root for the Jordan blocks of Jr which are not associ- 
ated with real negative eigenvalues, Kj and Kj are symmetric block diagonal 
matrices corresponding to the square roots of the real negative eigenvalues, 
and O is an arbitrary nonsingular symplectic matrix which commutes with 
diag(Jii, JjQ. 



Proof. Equation (|4.12l) is a special case of Equation (I4.11[) in Theorem 14.51 If 
Xj is a real square root of Jr then diag (Xj, — Xjj is a Hamiltonian square root 
of diag ( Jr, J^) and Hamiltonian structure is preserved under symplectic similarity 
transformations. There are 2 p+q real square roots of Jr which may be or not functions 
of Jr. 

For the second part, assume that Jr in (|2.17|l . 

C ni {ax,b\) 



Jr 



C np (dp, bp) 
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has only one real negative eigenvalue, say Xk < 0, k > p corresponding to the real 
Jordan block J nk . 

Let ±iM nk with M nk 6 R nx ™ be the two pure imaginary square roots of J„ fc , 
which are upper triangular Toeplitz matrices. See Corollary 14.11 and (|2.3p . Observe 
that (±iM nk ) 2 = —M 2 k = J nk . We will first construct a square root of diag ( J nk , Jj fc ) 
which is real and Hamiltonian. Let P„ fc be the reversal matrix of order which 
satisfies P 2 k = I (the anti-diagonal entries are all l's, the only nonzero entries). The 
matrices P rik M n . and M„ k P„ k are real symmetric and we have 



Mn k Pn k 



Pn k M rik 



—M 2 

■>ik 









PnkM nk Pn k _ 




L J n k \ 



Thus, 



Mn k Pn k 



~Pn k M n 



( and also 



-M nk Pn k 



P M 

r n k lll n k 



is a real Hamiltonian square root of diag ( J nk , J%.) which is not a function of 
diag(j„ fc , Jjj. 

If Xi = diag(F ni , . . . ,F np ) is a real square root of diag(C ni , . . . , C np ), 
X 2 = diag(L„ p+1 , . . . , £ ns ._i) is a real square root of diag(J„ p+1 , . . . , J ns ._i) and 
X 3 = diag(I/„ fc+1 , . . . ,L np+q ) is a real square root of diag( J nk+1 , . . . , J„ p+ J, then 



x 1 




x 2 




o 


MnkPnk 


x 3 










-xl 


~Pn k M nk 


o 




-xi 



X, Kj 
K 3 -Xj 



is a real Hamiltonian square root of diag (Jr, JjQ . Notice that there are 2 p+q different 
square roots with this form. Thus, for an arbitrary nonsingular symplectic matrix 
which commutes with diag (Jr, JjQ , 



y J (e) = *e 



Xj K$ 

k 3 -xj 



e- 1 *- 1 , j = i,. . . , 2 p+9 , 



is a Hamiltonian square root of W. 



If W has more than one real negative eigenvalue, the generalization is straight- 
forward. □ 
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5. Algorithms for computing square roots of a skew-Hamiltonian ma- 
trix. In this section we will present a structure-exploiting Schur method to compute 
a real skew-Hamiltonian or a real Hamiltonian square root of a real skew-Hamiltonian 
matrix W € R 2nx2n when W does not have real negative eigenvalues. 



5.1. Skew-Hamiltonian square roots. First we obtain the PVL decomposi- 
tion ofW £ R™ x " described in section [ 



U T WU = 



W x W 2 
O 



W4 = -W 2 , 



where U is symplectic-orthogonal and W\ is upper Hessenberg. The matrix U is 
constructed as a product of elementary symplectic-orthogonal matrices. These are 
the 2n x 2n Givens rotations matrices of the type 



cos ( 



sin 9 



sm t 



In—j 



1 < i < n, 



for some angle 9 6 [— ir/2, tt/2[, and the direct sum of two identical nxn Householder 
matrices 



Hj ® Hj(v, /3) 



In ~ (3VV 7 



where v is a vector of length n with its first j — 1 elements equal to zero. A simple 
combination of these transformations can be used to zero out entries in W to accom- 
plish the PVL form. See Algorithm 1 and Algorithm 5 in [31 pp. 4,10]. The product 
of the transformations used in the reductions is accumulated to form the matrix U. 

Then the standard QR algorithm is applied to W\ producing an orthogonal matrix 
Q and a quasi-upper triangular matrix Ni in real Schur form (|2.8p so that 



Wx = QNiQ T , 

and we attain the real skew-Hamiltonian Schur decomposition of W, 



T = WWU 



Ni N 2 ' 



No 



where U = U 



Q O 
O Q 



and N 2 = Q T W 2 Q. 
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This procedure takes only approximately a 20% of the computational cost the 
standard QR algorithm would require to compute the unstructured real Schur de- 
compositon of W [21 p. 10]. 

Let 



X Y 
X T 



Y = -Y 1 



be a skew-Hamiltonian square root of T. We can solve the equation Z 2 — T exploiting 
the structure. From 



'X Y ' 




'X Y ' 




~X 2 XY + YX T ~ 




'Ni N 2 










o (x T ) 2 




N*_ 



we have 



X 2 = N x 



(5.1) 



and 



XY + YX 1 



N, 



(5.2) 



Equation (|5.1|) can be solved using Higham's real Schur method (see Algorithm 1, 
Section 13.31 page [121 and it is not difficult to show that X inherits N\ 's quasi- upper 
triangular structure. Equation (15.2[) is a Lyapunov equation which can be solved 
efficiently since X is already in quasi-upper triangular real Schur form and Y is skew- 
symmetric. The techniques are the same as for the Sylvester equation. See [TJ Q31 
Chapter 16]. 



If the partitions of X = (Xij), Y = (Yij) and N 2 
block structure, 



(Nij) are conformal with N\ 



X 



X 



it 



A 12 
X22 



X\ m 
Xi m 

Xyy.™ 



Y 



Y 2 i 



-Y T 

Y 22 



1 ml 1 7r. 



-Y T ' 

ml 

-Y T 



Y 



No 



'N n 

N21 



22 



N, 



1 11 1 



N, 



hi 2 



N„ 



Y - V ' AT-- 



= -K 



1,.. 



11 ' 

rn. 
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then, from ()5.2[) . we have 



and 



k—i k—j 



XuYij + YijXjj — Nij — XikYkj — YfaX 



k=i+l 



k=j+i 



— Nij — ^ X ik Y k j — YkX 



k=i+l 



k=j+l 



(5.3) 



These equations may be solved successively for Y mm ,Y mjm -i, . . . ,Y m i, 
Y m -i,m-i,Y m -i tm -2, ■ • • , Vm-i,i) ■■■■> Y 22, Yii and Y n - We have to solve 



XuY^ +YijXfj — - XikYkj- Y lk Xj k + ^ Y ki Xj k , 



k=i+l 



k=j+l 



k=i+l 



i = m, m — 1, . . . , 1 



(5.4) 



Since Xa are of order 1 or 2, each system (|5.4j) is a linear system of order 1,2 or 4 
and is usually solved by Gaussian elimination with complete pivoting. The solution 
is unique because Xa and —Xj^ have no eigenvalues in common. See Section 

Algorithm 2 [Skew-Hamiltonian real Schur method] 

1. compute a real skew-Hamiltonian Schur decomposition of W , 



T = WWU 



'Nx N 2 
Nf 



2. use Algorithm 1 to compute a square root X of Ni, X 2 = Ni; 

3. solve the Sylvester equation XY + YX T — N 2 using (|5.4p and form 



X Y 
X T 



4. obtain the skew-Hamiltonian square root of N , X = UZU T . 

The cost of the real skew-Hamiltonian Schur method for W £ jj 2 ™* 2 " j s measured 
in flops as follows. The real skew-Hamiltonian Schur factorization of W costs about 
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3(2n) 3 flops [IT| |2]. The computation of X requires n 3 /6 flops, the computation of 
the skew-symmetric solution Y requires about n 3 flops [ITJ p. 368] and the formation 
of X = UZU T requires 3(2n) 3 /2 flops. The total cost is approximately 5(2n) 3 flops. 
Comparing with the overall cost of Algorithm 1, the unstructured real Schur method, 
which is about 17 x (2n) 3 flops, Algorithm 2 requires considerably less floating point 
operations. 



5.2. Hamiltonian square roots. Analogously, let Z be a Hamiltonian square 
root of T, 



Z = 



X Y 

-X T 



Y = Y 



(which is not a function of T). To solve the equation Z 2 = T, observe that, from 



X Y 




X Y 




X 2 XY - YX T ~ 




'Ax N 2 ~ 


-x T _ 




-x T _ 




o (x T ) 2 




Nl_ 



it follows 



and 



X = Ni 



XY-YX 1 =A 2 



(5.5) 
(5.6) 



Equation (|5.5p can be solved using Higham's real Schur method and Equation (JST 
is a singular Sylvester equation with infinitely many symmetric solutions. See [H 
Proposition 7]. Again, the structure can be exploited and we have to solve 



XuYij —YijXjj — Nij — ^ XikYkj + ^ Y lk Xj k - YkiXj k , 

k=i+l k=j + l k=i+l 

i = m, m — 1, . . . , 1 (5-7) 
j = i, i — 1, . . . , 1. 



The solution of the linear system (15. 7|) may not be unique but it always exists. 



Algorithm 3 [Hamiltonian real Schur method] 



1. compute a real skew-Hamiltonian Schur decomposition of W, 

Ai No' 
T=U J \VU= q N l 
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2. use Algorithm 1 to compute a square root X of N\, X 2 = AT X ; 

3. obtain one solution for the Sylvester equation XY — YX T = N2 using (|5.7p 
and form 



4. obtain the Hamiltonian square root of W, X = UZU T . 

6. Numerical examples. We implemented Algorithms 2 and 3 in Matlab 
7.5.0342 (R2007b) and used the Matrix Function Toolbox by Nick Higham available 
in Matlab Central website http : //www . mathworks . com/matlabcentral. To find the 
square root X in step 2 we used the function sqrtm_real of this toolbox and to solve the 
linear systems (|5.4|) in step 3 of Algorithm 2 we used the function sylvsol (the solution 
is always unique) . In step 3 of Algorithm 3 the linear systems (|5.7j) are solved using 
Matlab's function pinv which produces the solution with the smallest norm when the 
system has infinitely many solutions. 

Let X be an approximation to a square root of W and define the residual 



Then, we have X 2 = W +E and, as observed by Higham [T3j p. 418], the stability 
of an algorithm for computing a square root X of W corresponds to the residual E 
being small relative to W. Furthermore, for X computed with sqrtm_real, Higham 
gives the following error bound 



where || • || F is the Frobenius norm, c is a constant of order 1, n is the dimension of 
W and u is the roundoff unit. Therefore, the real Schur method is stable provided 
that the number 



Z = 



X 



Y 
X T 



E 



X 2 - W. 




a(X) 



11*1 



\\W\\ F 



is small. 



We expect our structure-preserving algorithms, Algorithm 2 (skew-Hamiltonian 
square root) and Algorithm 3 (Hamiltonian square root) to be as accurate as Algo- 
rithm 1 (real Schur method) which ignores the structure. The numerical examples 
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that follow illustrate that the three algorithms are all quite accurate when a(X) is 
small. 

Example 6.1. The skew-Hamiltonian matrix 



W = 



-A T 



A 








10~ 6 


1 










-1CT 6 





1 


icr 6 





A = 


-1 


-1 





10~ 6 


1 







-1(T 6 


-ifr 6 





1 










-l 


-l 






where e is the vector of all ones, has one complex conjugate eigenvalue pair and 3 
positive real eigenvalues (all with multiplicity 2). 

The relative residuals of both the skew-Hamiltonian and Hamiltonian square roots 
computed with Algorithm 2 and Algorithm 3 are 4 x 10~ 15 7 the same as for the square 
root delivered by Algorithm 1. 

Example 6.2. The eigenvalues of the skew-Hamiltonian matrix 



W = 



B 

A T 



A = 



B = 



' 




icr 6 







" 


icr 6 




























1(T 6 












icr 6 





" 


1 


2 


3" 






-l 





2 


3 






-2 


-2 





3 






CO 


-3 


-3 


0. 







are all very close to pure imaginary (four distint eigenvalues). 

The relative residuals of the square roots delivered by all the three methods are 
4 x 1(T 16 . 

If W has negative real eigenvalues there are no real square roots which are func- 
tions of W . However, all these algorithms can be applied and complex square roots 
will be obtained. In step 2 of Algorithms 2 and 3 a complex square root is computed 
and so we get a complex skew-Hamiltonian and a complex Hamiltonian square-root. 

Example 6.3. 

For random matrices A, B and C ( values drawn from a uniform distribution on 
the unit interval), the computed square roots of the skew-hamiltonian matrix of order 
2n = 50 (several cases) 



W = 



A 

C-C T 



B — B T 
A T 
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also have relative residuals of order at most 10 14 . 

7. Conclusions. Based in the real skew-Hamiltonian Jordan form, we gave a 
clear characterization of the square roots of a real skew-Hamiltonian matrix W. This 
includes those square roots which are functions of W and those which are not. Al- 
though the Jordan canonical form is the main theoretical tool in our analysis, it is not 
suitable for numerical computation. We have designed a method for the computation 
of square roots of such structured matrices. An important component of our method 
is the real Schur decomposition tailored for skew-Hamiltonian matrices, which has 
been used by others in solving problems different from ours. 

Our algorithm requires considerably less floating point operations (about 70% 
less) than the general real Schur method due to Higham. Furthermore, in numeri- 
cal experiments, our algorithm has produced results which are as accurate as those 
obtained with sqrtm_real. 
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